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Quantum Monte Carlo based on two-body density functional theory for fermionic 

many-body systems: application to 3 He 
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We construct a quantum Monte Carlo algorithm for interacting fermions using the two-body 
density as the fundamental quantity. The central idea is mapping the interacting fermionic system 
onto an auxiliary system of interacting bosons. The correction term is approximated using correlated 
wave-functions for the interacting system, resulting in an effective potential that represents the nodal 
surface. We calculate the properties of 3 He and find good agreement with experiment and with other 
theoretical work. In particular our results for the total energy agree well with other calculations 
where the same approximations were implemented but the standard quantum Monte Carlo algorithm 
was used. 

PACS numbers: 02.70.Ss,24.10.Cn,05.30.-d 



Density-functional theory 0, (DFT) and quantum 
Monte Carlo 0, (QMC) are generally thought of as 
two distinct approaches to the problem of interacting 
fermions. DFT is based on the Hohenberg-Kohn theo- 
rems 1] (HK) which state that the energy of an interact- 
ing fermion system in an external field can be written as 
a functional of the density, and that minimizing the en- 
ergy as a functional of the density gives the ground state 
energy (HK theorems) . Applications of DFT are usually 
based on the Kohn and Sham 0] method, where an aux- 
iliary noninteracting system is invoked. Minimization is 
achieved with respect to the orbitals of the auxiliary sys- 
tem. QMC also involves minimization of the energy. One 
way of minimizing, is to propagate a trial wavefunction in 
imaginary time [3|, so that it asymptotically approaches 
the ground state. 

In this letter we present a QMC method derived from 
an extension of DFT where the two-body density (n*- 2 )) 
is the fundamental quantity 0, E| • As in standard DFT 
the energy functional is universal but unknown, thus ap- 
proximations schemes are necessary. In the spirit of the 
Kohn-Sham ansatz we invoke an auxiliary system with 
identical as the system of interacting fermions un- 
der investigation, but instead of a non-interacting sys- 
tem, one of interacting bosons. As in the Kohn-Sham 
method, minimization is not performed with respect to 
the density, but with respect to the bosonic wavefunction 
via QMC Q| • This can be done, since the two-body den- 
sity can be written in terms of the bosonic wavefunction. 
In our method the sign-problem does not arise explic- 
itly, thus fixed-node 7] (FN) or released-node ]8] (RN) 
techniques are not needed. In the auxiliary fields Monte 
Carlo method 9] the need for FN or RN is also circum- 
vented, but the sign-problem still manifests in the phases 
of the auxiliary fields. 

The correction term necessitated by our ansatz is ob- 
tained approximately using correlated basis functions. 
The resulting approximation consists of a two-body and 
a three-body potential (effective nodal surface). The ap- 
pearance of the three-body potential (and density) in our 



energy functional is a result of our approximation scheme. 
In principle our approximate energy functional can still 
be written as a functional of n' 2 ', since according to the 
HK theorems the three-body density (as all other ob- 
servables) is a functional of rS 2 \ As far as the method 
developed here is concerned, the minimization itself is 
performed with respect to the bosonic wavefunction, thus 
higher-order potentials are easily handled. 

We apply our formalism to calculate the total energy, 
potential energy, and the structure factor of 3 He in a 
range of densities close to the equilibrium one (po = 
0.273A7 ct -3 , a = 2.556 A). Our model estimates the 
density that minimizes the total energy to be slightly 
less than the experimental result, as one would expect 
from the fact that we are not including back-flow effects. 
The calculated energies are in very good agreement with 
QMC results at the same level of approximation, and 
compare reasonably well with experiment. 

Given a system of interacting particles with potential w 
(including two-body and one-body potentials) and with 
two-body density (the diagonal elements of the two- 
body density matrix) HK can be extended as follows: 

• There is a one to one correspondence between w 
and n^ 2 \ 

• The ground state energy of the system can be ob- 
tained by minimizing E as a function of . 

The proof of these statements is an easy extension 
of the original work Q], and can be extended to 
A^-representable two-body densities using the Levy- 
constrained search [lfj. The energy (and all other ob- 



servables) can be written as a functional of rS 2 ^ as 
E[n {2) ] = 7> (2) ] + / drdr'w(r, r')n (2) (r, r'). 



(1) 



T[n( 2 )] is a universal functional of , that is, the kinetic 
energy can be determined by knowing only rS 2 " 1 . Whether 
the system is composed of bosons or fermions enters only 
in the form of T[n^]; the functional dependence of the 
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FIG. 1: Effective exchange pair potential. 

potential energy term on is the same for bosons or 
fermions. 

In the original DFT of HK, the universal functional in- 
cludes the kinetic energy and the pair interaction and is a 
function of the one-body density (-Fhk [«■])■ The analogue 
of -Fhk [n] in pair-density functional theory is simply the 
kinetic energy T[m 2 '], i.e. it does not include any of the 
potential energy terms. 

In order to obtain an applicable algorithm, we intro- 
duce an auxiliary system of interacting bosons, and add 
the required corrections. The issue of representability of 
a fermionic by a bosonic one shall be addressed in 
our extended study. Our starting equation is then 




where Te[n( 2 )] is the kinetic energy of a system of bosons 
with two-body density n^ 2 \ and 



AT[n (2) ] = T F [n (2) ] - T B [n (2) ], (3) 

where Tp[n^ 2 ^] is the kinetic energy of a system of 
fermions with two-body density nA 2 ' . 

The correction term in Eq. J2J is the difference be- 
tween the kinetic energies of two systems identical n/ 2 )'s. 
differing only in that one is a system of fermions, the 
other a system of bosons. We develop an approxima- 
tion to AT by constructing two trial wave-functions (one 
fermionic, one bosonic) with approximately equal n/ 2 )'s, 
and taking the difference of the kinetic energy expres- 
sions. The approximation presented here is valid for ho- 
mogeneous systems. 

For the fermionic one we take a wavefunction of the 
Jastrow-Slater form 

* F = D^F, (4) 



where 

F = \{f{n 3 ), (5) 

i<j 

and where (D^-) indicates a Slater determinant of 
plane waves between atoms of spin up(down), and /(r) is 
a correlation factor. In constructing a bosonic wavefunc- 
tion with the same two-body density we take the same 
correlation factor as in Eq. Q), but to account for the de- 
terminants we multiply by additional correlation factors 
between parallel spins, i.e. 

*b = FMF, (6) 

where 

^ m) =n/i a) M, a) 

where the product in Eq. 10 indicates a multiplication 

between pairs of parallel spins. The correlation factors 

t ( n 

fx {r) should be chosen in such a way that the two- 
body densities obtained from Eqs. and JHJ) are iden- 
tical. The correction term in this case can be explicitly 
obtained 

AT[n«] - To + J2 J dR\D^ D l \ 2 \7iF ■ V a T 

i 
i 

+E^v~ / ^ 2 (^iv 2 Fj^i),(8) 

i 

where R denotes all coordinates, m denotes the 
mass, and Nf,Nb are normalization integrals. To — 
3/5fcp/(2m), which is the kinetic energy of the homo- 
geneous non-interacting system (fcp is the Fermi wave 
vector). If the correlation factor f x is chosen such that 
the two-body and three-body densities obtained from the 
determinants in Eq. 10} are identical to those obtained 
from the correlation factors f x in Eq. © then the second 
and third terms in Eq. JHJl cancel resulting in 

AT[„( 2 >] = To + ]T J rfRF 2 (TjFi V 2 TjFi), 

(9) 

In obtaining a first approximation to the correlation 
term f x in the case of a homogeneous system we can 
make use of the radial distribution function of the non- 
interacting fermion gas, given by 

g x (r) = 1 — < „ (sin/cpr — kprcoskFr) 1 (10) 
^ Kpr J 

We obtain f x from an inverse hypernetted chain equation 
(we also tried using f x = g x , and obtained very similar 
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FIG. 2: Total energy for 3 He as a function of density. The 
solid curve is a fit to experimental results (Ref. fTgh . the curve 
connecting our calculated points is a fifth order polynomial fit. 
The triangle is the result from a fixed node DMC calculation 
without backflow 



results). While fx does not guarantee that the second 
and third terms in Eq. (JHJ will cancel, in this work we 
assume Eq. as the form for our approximation. Thus 
the correction term resulting from Eq. iJJJJ in our scheme 
is the sum of a two-body and a three-body interaction, 



1 |d 2 ln/,(r) , 2d\nf x {r) 
2m 
1 



dr 2 



Or 



w^(r 12 ,r 13 ) = — (Viln/^ria) • Viln/ X (n 3 )), (11) 
Ira 

and a constant term To. The expressions are the same 
for both /J, f*, however in general /J ^ /J. Since in 
this work we will deal with the unpolarized case, where 
fx — fx fr° m n ow on we will drop the arrows. 

Going beyond the weak coupling approximation would 
lead to correction terms including the correlation factor 
/. In principle / is a functional of the two-body density 
of the system, thus a self-consistent algorithm would be 
necessary. Possibly this can be avoided by using known 
correlation factors for a given system under investigation 
(or obtaining one from solving the Eulcr equation 11]). 
Potentially, better approximations can also be obtained 
for AT if better wavefunctions are chosen in Eqs. 10} and 
(|rj|) . In our scheme three-body correlations and Feynman- 
Cohen back-flow |l2( have not been considered. 

As a test of the quality of f x we have performed 
a Monte Carlo simulation at the experimental density 
(po = 0.273-/V/(7~ 3 ), and have found that the radial dis- 
tribution function obtained is in good agreement with 
that of the noninteracting Fermi gas. Obviously, in this 
case the energies between the fermionic system and the 
auxiliary bosonic one correspond. In Fig. ^ the effec- 
tive pair potential that incorporates the nodal surface 

(2) 

(-wi ; ) is shown at po- The potential is an effective way that the nedlect of back- flow corrections leads to an over- 



TABLE I: Energy quantities per particle as a function of 
density. All data are in Kelvin. 



of including the nodal structure, it is repulsive at short 
distances, and displays alternating valleys and barriers of 
decreasing magnitude. 

We applied the above procedure to a system of un- 
polarized 3 He atoms at zero temperature. We have 
used the HFD-HE2 interaction potential due to Aziz et 
al 13]. Between particles with parallel spin the poten- 
tial interaction modified by the additive terms given in 
Eq. ©. To reduce the variance in DMC we have used 
a guiding function of the form given in Eq. © with 
f(r) = exp - b 5 /2r 5 (b = 1.15a, where a = 2.556 A). 
The parameter b was optimized by a variational Monte 
Carlo calculation. For other calculations on the same 
system see Refs. [H Q Ql El 03 

We perform a series of calculations using the standard 
bosonic DMC0 algorithm. Our cell included 108 parti- 
cles in all cases, we used an imaginary time step of 50 a.u. 
We collected averages over 100,000 steps. Total energies 
are estimated in the standard way, coordinate dependent 
observables were estimated using pure estimators 0. 
The non-coordinate dependent part of the fermionic ki- 
netic energy was obtained by subtracting from the total 
energy the potential and the correction term in Eq. I|ll(l . 

In Fig. |2Jwe compare our calculated total energies with 
experimental results. The thick solid curve is a fit to ex- 
perimental results fl9| . The minimum density obtained 
by us is in close agreement with the experimental result 
(exp:0.273iVCT- 3 ,calc:0.244Arcr- 3 ). The structure factor, 
shown in Fig. [3] also compares well with experiment. 
The experimental [2(| minimum energy is —2.473 K, our 
calculated energy at that density is —2.147(7) K. Our en- 
ergy at the calculated minimum density (from the fitted 
curve) is —2.220(6) K. It is important to note that our 
results are in good agreement with previous calculations 
that do not include the back-flow correction (in Ref. ^3 
a QMC calculation without back-flow is reported result- 
ing in —2.128(15) K for the total energy per particle). 
We are not aware of fixed-node DMC calculations with- 
out back-flow for the full density range presented here, 
but variational Monte Carlo |TlL fl5| calculations indicate 
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FIG. 3: Comparison of calculated (solid triangles) and ex- 
perimental (circles: Ref. |24| and empty triangles: Ref. I25T) 
structure factor of 3 He at equilibrium density. 

estimation of the energy by a few tenths of a Kelvin (see 
Figure 1. in Ref. Ilojh Our calculated energies differ more 
from the experimental ones at higher densities. This is 
not surprising, since it is known that the back-flow ap- 
proximation is more crucial at higher densities. |2lj 

In Table[l]we present values for the total, potential, ki- 
netic energies, and for the coordinate dependent part of 
the correction term (AT — To). Our value for the poten- 
tial energy at po (—14.03(3) K) also compares reasonably 
well with other theoretical results ((V) = -14.84(10) K 
[lfij ]. The coordinate dependent part of AT gives only a 
small correction compared to the value of the kinetic en- 
ergy itself (the bosonic kinetic energy plus To is already 
a reasonable approximation to the kinetic energy). Thus 
using an auxiliary bosonic system is a promising scheme 
for developing approximations. 

We have demonstrated that using an algorithm con- 
structed from the pair-density a good description of an 
interacting fermionic system can be obtained. Our al- 
gorithm is arrived at by invoking an auxiliary system of 
bosons, therefore the calculation itself can be performed 
by a bosonic DMC algorithm. While it is clear that fur- 
ther work needs to be done to obtain a better approxima- 
tion, the fact that we have obtained quantitative results 
for the observables calculated demonstrates that this av- 
enue is worth pursuing. Our future directions include the 
developing and testing of more sophisticated approxima- 
tions for the kinetic energy correction term used here, 
such as implementing back-flow, three-body correlation, 



and inverti ng t he fermionic hypernetted-chain approxi- 
mation J22 L |23| . Comparison of our method to related 
ones is also of interest. 

This research was supported by MIUR-2001/025/498 
and by SISSA. We benefitted greatly from discussions 
with Professor K. E. Schmidt. 
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